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ABSTRACT 

We study the spherical Mexican hat wavelet (SMHW) as a detector of primordial 
non-Gaussianity of the local type on the Cosmic Microwave Background (CMB) 
anisotropies. For this purpose we define third order statistics based on the wavelet 
coefficient maps and the original map. We find the dependence of these statistics in 
terms of the non-linear coupling parameter /„/ and the bispectrum of this type of 
non-Gaussianity. We compare the analytical values for these statistics with the re- 
sults obtained with non-Gaussian simulations for an ideal full-sky CMB experiment 
without noise. We study the power of this method to detect f n i, i. e. the variance of 
this parameter a 2 (f n i), and compare it with the variance obtained from the primary 
bispectrum for the same experiment. Finally we apply our wavelet based estimator on 
WMAP-like maps with incomplete sky and inhomogeneous noise and compare with 
the optimal bispectrum estimator. The results show that the wavelet cubic statis- 
tics are as efficient as the bispectrum as optimal detectors of this type of primordial 
non-Gaussianity. 
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1 INTRODUCTION 

The primordial perturbations generated during the inflation- 
ary period are imprinted in the radiation and matter distri- 
bution. The study of the CMB anisotropies has become an 
important source of information to understand the physics 
of the very early universe. Thus for example the search of 
primordial non-Gaussianities on the CMB anisotropies, the 
measurement of the tilt and running of the index of the 
power spectrum A^(fc) = Ak" 3 ^ 1 and the search of pri- 
mordial gravitational waves have become part of a set of 
observables that are being used to select among many dif- 
ferent models for the inflationary paradigm. The simplest 
models of inflation as the standard, single-field, slow roll in- 
flation JGuthll98ll;lAlbrecht fc Steinhardtl 19821 ; |LincUll982t 



on the Wilkinson Microwave Anisotropy Probe WMABJ 
data, there are several reports of anomalies present 
in these data. We can mention the as ymmetry be- 
tween the two e c liptic hemispheres (lEriksen et al.l 
20041; lHansen et alj |2004 [Eriksen et al.l 120051. 120071; 



Hoftuft et al.l 120091 : iPietrobon et all |2010| ; [Vielva fc Sairj 
2010 bJ), anomalous quadrupole-octupole alignment 



JCopi et al.ll2004l; Ide Oliveira-Costa et alJl20bH; ICopi et al 



20061 ; iGruppuso fc Buriganal 12009] ; iFrommert fc E n Him 
20ldh a non-Gaussian cold spot (I Vielva et al 



. 12004 ; 
Mukheriee fc Wand 12004 ICruz et all 120051 . l200d . l2007al lb. 



20081 ; IVielyaet al.N2010af), unexpected alignme nt of CMB 



structures |Wiaux et alj|200q ; I Vielva et al,l2007|) and an un 



1 19831 ; iMukhanov et al.l |l992j) predict that the anisotropies 
are compatible with a nearly Gaussian random field. 

There are two main procedures for the Gaussianity anal- 
yses of a map of CMB anisotropies. One can perform blind 
tests searching for deviations with respect to the null hy- 
pothesis (the random field is Gaussian). The second option 
is to consider different physical scenarios and to look for 
their imprints on the anisotropies. 

Among the many different blind tests performed 



expe c ted low value of the CMB variance i Monteserfn et al.l 
120081 ; ICruz et alj|20ich . 

Regarding targeted tests, one has to think about 
the possible physical mechanisms that lead to non- 
Gaussianities and search for their possible signatures. His- 
torically motivated, many inflationary models that gener- 
ate non-Gaussianity can be parametrised by the local non- 
linear coupling parameter /„;, which is i ntroduced through 
the primordial gravitational potential |Verde et al.l I2OOO1 ; 
iKomatsu fc Spergelll200ll ; lBartolo et alj|2004h ~ 
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$(x) = $ £ (x) + f nl {<f>l (x) - («!>l(x)) 2 }. 



(1) 



Significant non-Gaussianity of the local form c an be gen- 
erate d for example in the curvato n model (ILvth et al] 
|2003| >. multi-field inflationary models ilKomatsu et al] 20051. 
the inhomogeneous reheating scenario i Dvali et al] 12004 ; 
Bart olo et al] 12004), models with low reh e ating efficiency 
l|Mukhanov et al] ll99J ISalopek fc Bond! ll99oh . models 
based on hybrid inflation (|Linl |2009j), etc. This kind of 
non-Gaussianity is characterised by correlations among 
modes k in the Fourier space with very different amplitude 
IjCreminelli et al.ll2006l . 120071 1. This can be seen for example 
through the shape function F(ki,k2,kz) of the local non- 
Gaussianity 



F(k 1 ,k 2 ,k 3 )=Af nl 



+ 



f\i-\ Kin fCi A/q 



+ 






(2) 



where A is a normalization constant (see for example the 
plot of the shape function of the local distribution by 
iBabich et al] 12004 1. The most significant contributions of 
this kind of non-Gaussianity arise for the cases with k\ « 
&2 ~ &3 and permutations among the three modes fci, 
ki and k%. Other shapes of inflationary models that pro- 
duce their particular kind of non-Gaussianit y are for exam- 
ple the equilateral and the orthogonal shape (jSenatore et al] 
l2010l : iKomatsu et al.ll2010l ). The most significant contribu- 
tion to the non-Gaussianity is located in specific ranges of 
the Fourier space k\ ~ k-2 ~ k% for the equilateral shape 
whereas the orthogonal shape is nearly orthogonal to the 
two previous forms. 

In this paper we focus on the non-Gaussianity of 
local type. The canonical inflationary model predicts 
fni ~ 10~ 2 whereas other models predict larger amounts 
of non-Gaussianity (jBartolo et al] 120041) . An eventual 
detection of a deviation from Gaussianity will rule out 
many inflationary models from the present available list. 
Many studies have been performed to constrain the local 
f n i on different data set s such as the Cosmi c Background 



Explorer (COBE) data (Komatsu ct al. 2002; Cayon et al 



2003a|), the MAXIMA data (|Santos et a l. 2003; Cayo n et al 



the B OOMERang data (|De Troia eTal . 2007; Natoli et all 



2003bl l. the Very Small Array data iSmith et al] l2004h. 



|2009| 1 and the Archeops data (|Curto et al] 120071 . l200Sft . 
Significant improvements on the (S/N) ratio of the local 
/„; parameter have been achieved with the WMAP data 
using different estimators. For most of the works, /„; is 
positive with a significance between la and 2a. See for 
example the results obtaine d with differe n t bi spcctrum- 
based e stimators (Komats u et al] 120031; [Spergel et al] 



20071; 



200' 



2010) 



_ . , Creminelli et al. 2006, 2007; Yadav & Wandclt 2008; 
Komatsu et al.l 120091: ISmith et ail 120091; [Eisner fc WandeltJ 
iBucher et al.ll201(il ; IKomatsu et alj|20icl: ISmidt et al] 
wit h the SMHW (ICurto et al] l2009allbf) , with 
needl e ts dPietrobon et al] 120091 ; iRudiord etail 120091, 
l20ld; ICabella et al] l201Ch . with the HEALPix wavelet 



I Casapon sa et al] 2010J). with the M i nkowski functionals 



( Hikage et al] |2006| ; iGott et al] 120071 : iHikage et al] |200S|; 
Matsubarall20ld 1. the N-PDF distribution (IVielva fc Sanzl 



2009), the skewness of the power spectrum (Sm idt et al 
20091 ), etc. Other works use the dis tribution o f mat- 
ter on large scales (see for example I Dalai et al 
iMatarrese fc Verdei [20081 ; ISlosar et al] 120081 ; ISeliai 



2008; 



2009 



iDesiaccmes fc Seliakll2010l ; IXia et al]|2010l 1 to constrain the 
local f n l. 

This work is a study of the effici ency of the third or- 
der wavelet-based estimators (jCurto et al] l2009al lbr) to de- 
tect primordial non-Gaussianity of the local type. We find 
the dependence of the third order estimators in terms of 
the local bispectrum and f n i . We compare the power of this 
method to detect f„i with the optimal bispectrum estima- 
tor for an ideal full-sky and noiseless experiment and the 
same comparison for an experiment with WMAP 5-year and 
WMAP 7-year beam, noise and sky cut properties. Our re- 
sults indicate that the wavelet estimators are as efficient as 
the optimal estimators based on the bispectrum. 

The article is organised as follows. Section [2] presents 
the estimators based on wavelets, their dependence on the 
angular bispectrum, their analytical covariance matrix and 
the /„; Fisher matrix for wavelet and bispectrum estimators. 
In Section [3] we compare the a(f n i) values obtained with 
the wavelet estimator and the bispectrum estimator and the 
conclusions are presented in Section [4] 



2 THE THIRD ORDER STATISTICS 

2.1 Expected values of the wavelet estimator in 
terms of the bispectrum 

The third order stati stics of this analysis ar e base d 



1998 



2007); 



on the SMHW. See lAntoine fc Vanderghevnstl 
Martinez-Gonz alez et al] ( 20021 ) ; IVielva 

Martinez-Gonzalezl ( 20081 ) for detailed information about 
the wavelets and a list of applications to the CMB 
anisotropics. Given a function /(n) defined at a position n 
on the sphere and a continuous wavelet family on that space 
^(n; b, R), we define the continuous wavelet transform as 



(i?;b) = / dn/(n)*(n;b,i?) 



(3) 



where b is the position on the sky at which the wavelet 
coefficient is evaluated and R is the scale of the wavelet. 

Considering a set of different angular scales {Ri} we de- 
fine a third order st atistic depending on three scales {i, j, k} 
jCurto et alj|2009bl ) 



Qijk 



Ail aiOjau 



dnw(Ri, n)w(Rj, n)w(Rk, n) 



(4) 



where Oi is the dispersion of the wavelet coefficient map 
w(Ri,n). In the particular case of Ro = 0, w(Ro,n) = /(n). 
Using the properties of the wavelet, we have 



■](Ri,n) = 2 atmU)t(.Ri)Yt m (n) 



and 



E^< 



21+1 2, R , 



(5) 



(6) 



where uje(R) is the window function of the wavelet at a 
scale R and it is given by the h armonic transform of the 
mother wavelet of th e SMHW (|Martinez-Gonzalez et al] 
l2002l ; ISanz et alj|200d l. The convolution with the wavelet is 
equivalent to filter the maps with a window function U)t(R) 
which depends on the scale. In Fig. [T] we plot the wavelet 
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Figure 1. The window function for the SMHW at different an- 
gular scales. Note that the wavelet filters high multipoles for low 
angular scales and viceversa. 



window function for several angular scales. For small scales, 
the wavelet niters low multipoles and viceversa. Therefore 
it is important to select a set of angular scales that ranges 
all the interesting multipoles. 

2.2 The statistics and the primordial 
non-Gaussianity 

Considering Eqs. [4] and [5] the third order moments can be 
written as 



Qijk — 



4-7T UiajCTk 
%mi "£ 2 m 2 ^^3 m 3 



1 £ 

V. ^1,^2^3, m l 1^2,7713 

xu(, {Ri)ui 2 {Rj)uJt 3 (R k ) 
*i 4 4 \ ( *i £2 4 \ 



77li m.2 TO3 







(2& + 1)(24 + 1) (24 + 1) 



4tt 



(7) 



where matrix is the Wigner 3j symbol and we have used the 
Gaunt integral (jKomatsu fc Spergelll200 u ) 



d nYe imi (n)Ye 2m2 (n)Yi 3Tn3 (n) 
£2 £3 \ ( £1 £2 £3 



mi rri2 rns 







(2*i + l) (2^ + 1) (24 + 1) 



47T 



(8) 



The mean value of the third order statistic gyj can be 
written in terms of the re duced bispectrum as defined by 
iKomatsu fc Spergell (|200ll ): 



(qijk) 



4-K OiOjOh 
2 , Ul l (Ri)Vt 2 (Rj)Ui 3 {Rk) I I 1 t 2 t 3 bl 1 l 2 l 3 , 



(9) 



where Ieii 2 e 3 is defined as 







4 4 





(2*i + l) (2*2 + 1) (24 + 1) 



4tt 



(10) 



Now assuming a primordial gravitational potential of the 
form given by Eq. [T|it is possible to derive its correspondin g 
bispectrum in terms of f n i fsee IKomatsu fc Spergell 1200 11 ). 
The expected value of the third order moments is propor- 
tional to fni 



(««*)/»! = Qijfe x f n i, 
where 



(11) 



*ijk — 



47T (JiOjOk 

'{ ^2 w «i(- Ri )^ 2 (-Rj)wi 3 ( J Rfe)i'| 1 £ 2 i 3 &?^ 3 |. 



(12) 



The pixel properties 
placing Ce by Ci 



are taken 
2 

in Eq 






(pix) (pix) (pix) 
W B . ^tf„ 0J/>„ 



into account by 
El and 6f^ 3 

(pix) 



by 



in Eq. 1121 where lo\ is the 

pixel window fun ction for the HEALPix pixelization 
l|G6rski et alj|2005h . 

We have evaluated the ayjb statistics analytically us- 
ing Eq.[l2]as well as from non-Gauss ian simulations for t he 
same set of 12 angular scales used in lCurto et al.l i|2009bl l 2 l. 
We have used a full-sky ideal experiment without noise and 
a characteristic angular resolution of 6.9 arcmin (HEALPix 
N S ide = 512). We need the expected values of the local 
primordial bispectrum to evaluate analytically the statis- 
tics aijk- We have computed the primordial bispectrum up 
to £ m ax = 1535 using the gTfa st|j code based on CMB- 
Fast (|Seliak fc Zaldarriagal Il996f ) to evaluate the transfer 
function. The cosmological parameters for this analysis are 
Q. cdm = 0.25, n b = 0.05, n A = 0.70, r = 0.09, h = 0.73 
and a scale invariant spectral index n — 1 for the power 
spectrum P(k). We have used a set of 300 non-Gaussian 
simulations generated with the same c osmological p a rame- 
ters f ollowing the algorithm defined in lLiguori et al.l (|2003l . 
120071 ). The mean value of the aijk statistics of these simu- 
lations and its error bars are plotted in Fig. [5] These values 
are compared with the theoretical value obtained from Eq. 
1121 We can see that there is a good agreement between the 
simulations and the analytic model. There is a slight dis- 
crepancy at negative values but we do not consider it very 
significant due to the correlations among the aijk quanti- 
ties. In any case, these values correspond to CHjk involving 
large scales where cosmic variance is more important and 
a larger number of simulations would be needed to achieve 
convergence. There is also a slight deviation at small angular 
scales (corresponding to values of aijk ~ 0.5 x 10 -4 ). This 
may be related to numerical uncertainties in the evaluation 
of the integrals that lead to the bispectrum &J?"™ at high 
multipoles *. However, we have seen that these small devia- 
tions do not introduce significant differences (< 1%) in the 
estimation of er(fni) when we use either the analytic or the 
simulated oyh. 



Uxt2la 



2 The angular scales used in lCurto et al.l J2009n ) are: 6.9, 10.6, 
16.3, 24.9, 38.3, 58.7, 90.1, 138.3, 212.3, 325.8 and 500 arcmin. 
The unconvolved map was also included in the analysis. 

3 http://gyudon.as.utexas.edu/~komatsu/CRL/nongaussianity/ 
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Figure 2. The a^k statistics computed analytically and with 
nsim = 300 non-Gaussian simulations. The error bars correspond 
to the dispersion tj(cnjk )/V re sim obtained with the simulations. 
The ctjjfc are sorted such as the largest values correspond to 
pixel-dominated scales while the smallest values correspond to 
the largest scales. 
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Figure 3. The analytical covariancc matrix computed through 
Eg. 1 161 compared to the one obtained from Gaussian simulations 
for the same power spectrum Ci. 



2.3 The covariance of the third order statistics 

The covariance matrix of the third order moments can be 
computed in the Gaussian limit using the properties of the 
covariance for the bispectrum. Considering the Gaussian 
case we have 



Cijk,rst — (qijkQrst) ~ \Qijk) \Qrst) — (qijklrst)- 

From the definition of qijk 

1 1 1 



(13) 



(qijhQrat) — 



X 



dnidn 2 



(47r) 2 OiOjOk O r O a O t 

x (w(Ri, m)w(Rj, ni)w(Rk,n±) 
xw(R r , n 2 )w(R 3 , n 2 )w(R t , n 2 )). 



(14) 



Using Wick's theorem and the properties of Gaussian distri - 
butions (in a similar manner as in Eq. 13 of iHeavensi 1 19981 ) 
we have 

(w(Ri, n 1 )w(R j ,n 1 )w(R k ,n 1 ) 

xw(R r ,n 2 )w(Rs,n 2 )w(R t ,n 2 )} = 

(w(Ri,nx)w(Rj,ni)) 

x {w(R k , m)w(Rr, n 2 )> 



x (w(R s , n 2 )w(R t , n 2 )) + 
+ permutations (total 15 terms) 



(15) 



From these terms, there are only 6 that do not vanish in Eq. 
1141 They are those which in volve the two co ordinates ni 
and n 2 in the same average (|Srednickilll993l ). Taking this 
into account and the pr operties of the two point correlation 
functions (see Eq. 14 of lHeavens!ll99ct) we have 



(qijkQrst) — 



(47r) 2 GiGjOk o- r a s a t 



7 t Ihi 2 i 3 ChCi 2 Ci 3 

x {u]i ± (R l )uj h (Rr)oJi 2 (Rj)ui 2 (Rs)^i 3 (Rk)u>i 3 (R t ) 
+ui 1 (Ri)ui x (R r )uji 2 (Rj)u)i 2 (R t )u)i 3 (Rk)toi 3 (R s ) 
+ui 1 (Ri)oJi 1 (R a )ui 2 (Rj)ui 2 (Rr)ui 3 (Rk)ui 3 (R t ) 
+ui 1 (Ri)ui x (R B )ui 2 (Rj)loi 2 (R t )uji 3 (Rk)toi 3 (R r ) 
+uj Li (Ri)ui x (Rt)coi 2 (Rj)uji 2 (R r )uJi 3 (Rk)toi 3 (R s ) 

+0J tl (R t )uj h (Rt)(Jl 3 (Rj)0Jl 2 (Rs)0Jl 3 (Rk)0Jl 3 (Rr)} 

(16) 

The pixel properties are taken into account here by replac- 
ing Ci by CiU£ u^ . We have compared the analytical 
covariance matrix obtained through Eq. [16] with the covari- 
ance matrix obtained with 10 4 Gaussian simulations for the 
same cosmological parameters defined in Subsect. 12.21 The 
covariance matrix elements djk,rst, are compared by pairs 
in Fig. [3] obtaining a very good agreement. 



2.4 /„( Fisher matrix of the third order moments 

We discuss the detectability of primary non-Gaussianity 
with the third order moments. Assuming that the third or- 
der statistics are Gaussian-like, we can use the Gaussian 
likelihood to constrain the / n j parameter 



L(fni) = C e 



-X^(/„()/2 



(17) 



where Co is a constant and X 2 (fni) is given by 



X 2 (fnl) = ^ (liJk-(liJk)f nl )C i / ktrgt (q° b s S t -(qrst)f nl )-(l8) 
ijk,rst 

C~l . is the inverse of the covariance matrix of the third 

zjk,rst 

order statistics that we have computed analytically in the 
previous section, q°^ k axe the third order statistics obtained 
from the data and (qijk)f nl ar e the expected values of the 
third order statistics for a given model with /„; . As we have 
seen in Subsect. 12. 21 (qijk)f nl = fnlOCijk, with a^k a constant 
independent of f n i ■ Using this on Eq. 1181 we have 

X 2 (fnl) = ^ fay* - a »jft/n;)C l ^fe,r S t( l ?rL S ~ Oirstfnl)- (19) 

ijk,rst 

The variance of the f„i parameter can be computed using 
the Fisher matrix, which leads to: 



a (fm) = 



-1 



B2logL{f nl ) 



<>r 



i 9 2 x 2 (f n i) 

1 



l~iijk,rst ai 3 k( - y ijk,rst a rst 



(20) 
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2.5 Principal Component Analysis of the third 
order moments 

One of the most significant advantages of the wavelet-based 
analysis for this type of non-Gaussianity is that we are able 
to reduce the non-Gaussian information present in the bis- 
pectrum b( 1 i 2 e 3 (about 10 7 elements) to a set of several hun- 
dreds of qijk statistics. However, these quantities are corre- 
lated and in certain conditions their covariance matrix may 
be ill-conditioned. There are different approaches to deal 
with these matrices. We will apply a Principal Component 
Analysis (PCA) in all the tests where the wavelet coefficient 
covariance matrix is involved in order to keep the errors re- 
lated to the covariance matrix under a certain threshold. 

The x 2 statistic of a random variable x of dimension n 
can be written 




1000 



Figure 4. o"(/ n ;) for different £ m ax using the Fisher matrix of 
the bispectrum. 



- 2__, X iCij J 



(21) 



»>j=i 



where C is the covariance matrix of this variable. For any 
positive definite covariance matrix, it is possible to find a lin- 
ear transformation of the x vector where the corresponding 
covariance matrix is the identity 



2 V^ 2 

x =2^Vi 



(22) 



where yi = y^"_ (D ' R^jjXj and D and R are the eigen- 
value and eigenvector matrices of the covariance matrix, 
C — RDR f . In particular we sort the eigenvalues in descend- 
ing order Di > Di+%. We can define a partial Xm statistic 



2 _ V"^ 2 

Xm / Vi 



(23) 



such that only the yi quantities related to the largest m 
eigenvalues are considered. The Fisher matrix for the f n i 
parameter using Xm is defined as 



°~m{fnl) 



1 dxLUnl) 

2 dJT, 



(24) 



As the third order moments are linearly proportional to f n i , 
x j = li — fnldj and defining /3i = J^™ (-D -1 / 2 ^^^, tne 
Fisher matrix for the /„,; parameter is 



Cm(/»l) = 



TZ^f 



(25) 



The unbiased bispectrum-based minimum variance estima- 
tor for the incomplete sky limit and inhomogeneous noise is 



tor tor the inco mplete si 
l|Creminelli et alj|2006h 



* nl ~ at Z^ 



&£\m,\ 6^2 m 2 ^3 m 3 



),,-, 



xC, 



.c; 



.c: 






(27) 



In particular, the signal-to-noise ratio for the /„; parameter, 
using a Fisher anal ysis for a full-sky and homogeneous noise 
experiment is IjKomatsu fc Spergelll200ll ) 



fnl 



jjrnodel rtmodel 

&t 1 e 2 t 3 - ti i 1 t 2 e 3 



(28) 



2^e 1 ^e 2 ^e 3 *i*2*3 

where o^ t t is the variance of the bispectrum. 

Given a map with N P i X pixels, the number of operations 
to obtain the local / n ; through the KSW algorithm based 
on the bispectrum is ~ 100 A^ 2 (JKomatsu et al.ll2005l ). Our 

algorithm based on the SMHW requires ~ n sca iN '( op- 
erations, where n sca i is the considered number of scales 
(n sca i ~ 15 for a WMAP-like experiment). Therefore, the 
wavelet algorithm is about an order of magnitude faster than 
the bispectrum algorithm for this kind of non-Gaussianity. 



2.6 Bispectrum estimator and error bars 

The bispectrum-based estimators are the most widely ap- 
plied tools for detecting primordial non-Gaussianity. Con- 
sidering the angle averaged bispectrum Be 1 t 2 £ 3 as defined in 
iKomatsu fc Spergell l|200ll ). it can be shown that the unbi- 
ased bispectrum-based minimum variance estimator for the 
full-sky linrit and homogeneous noise is (|Creminelli et al.l 
l2006h 



f = I V ( ix e ' 2 iz \ BT ^ 

Jnl AT Z^ I mi m 2 m 3 J C ei C'e 2 Ce 3 



(26) 



3 COMPARISON WITH THE OPTIMAL 
BISPECTRUM ESTIMATOR 

In this Section we study the values of cr(fni) obtained with 
the wavelet and bispectrum estimators for a full-sky and 
homogeneous noise experiment and for WMAP-like 5-year 
and 7-year V+W combined maps. 



3.1 Estimated error bars for / n j with an ideal 
experiment 

We consider the local primordial bispectrum of a noiseless 
experiment with an angular resolution of 6.9 arc minutes 
(imax — 1535), and a cosmological model characterised by 
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Figure 5. The window function of the SMHW for the list of selected angular scales for this case. We have convolved each wavelet window 
function with the pixel window function at N sia - e = 512. Note that the scales of 2.9 and 4.5 arc minutes are smaller than the pixel size 
(6.9 arc minutes). However they provide additional information because of the shape of the wavelet. 



Sl cdm = 0.25, Q b = 0.05, Ov = 0.70, r = 0.09, h = 0.73 and 
n — 1. The variance of /„; is computed through the Fisher 
matrix of the bispectrum (see Eq. I28f> . We have used the 
gTfast code to estimate the radiation transfer function and 
evaluate the expected values of the local bispectrum for thi s 
experiment following Eq. (19) o nKomatsu fc Spergell IJ2001T) . 
We have computed <t(/„j) for different (z ^ (.max (see Fig. 

S). 

Note that as we include higher multipoles (i.e. smaller 
scales) <r(f n i) decreases obtaining a limit of <r(f n i) = 5.4 for 

*-ma.x — lOoO. 

For the wavelet estimator, we consider the ideal exper- 
iment and the cosmological parameters defined previously 
for the primordial bispectrum estimators. To maximise the 
non-Gaussian signal we take a wide interval of different 
angular scales. We select 21 angular scales (including the 
unconvolved map) from 2.9 arc minutes to 167.07 degrees 
logarithmically spacecy. See Fig. \5\ for the list of consid- 
ered scales and the wavelet window function correspond- 
ing to these scales. We compute the third order quantities 
ctijk given by Eq. [12] and the covariance matrix of the qijk 
statistics given by Eq.[TS]for these angular scales. There are 



4 The two smallest scales peak at very small angular scales, cor- 
responding to sub-pixel structures at the selected resolution of 
N s ide = 512. Although the localization properties of the wavelet 
may be affected for those scales for the considered resolution, we 
have checked that this effect does not affect the /„; estimation. 
To test this point, we have compared the values of <r(/ n ;) using 
the same set of 21 angular scales in the current case of pixel re- 
sultion N aia - e = 512 and £ m ax = 1535 and in another case with 
N s ide = 1024 and (max = 1535. The differences in <x(/ n ;) are 
< 1% between both cases. 



(21+3-l)!/[3!(21-l)!] = 1771 different third order statistics 
for 21 angular scales. The variance of f n i given by Eq. [20] 
and the variance estimated with Monte Carlo simulations 
require the inverse of the covariance matrix C'ijk,rst- This 
covariance matrix has a large condition number, defined as 
the ratio of the maximum and minimum eigenvalues. The 
eigenvalues of this matrix are plotted on the left panel of 
Fig. [6] This implies that the computation of its inverse is 
an ill-conditioned problem. Therefore the numerical errors 
present due to the limited precision of our computers (of 
the order of 2 1 " 53 ~ 10" 16 ) can affect the value of a(f„i). 
In order to take into account the effect of this source of 
errors, we compute a(f n i) for different subsets of eigenval- 
ues of the C matrix using the PCA described in Subsect. 
12.51 This is plotted on the right panel of Fig. [F5] using the 
Fisher matrix for /„; given by Eqs, 1241 and 1251 and the PCA 
method. We as well plot the value of <r(f n l) obtained with 
the Fisher matrix of the bispectrum given by Eq. [25]for the 
same experiment and (.max = 1535. In particular, imposing 
a limit on the ratio of the minimum and maximum eigen- 
values of the covariance matrix Di/D max ^ 10~ 12 , we have 
cr(fni) =5.4. Similar results are obtained for a(f n i) when it 
is estimated with 10 3 Monte Carlo simulations. In particu- 
lar, we have obtained a Fisher (f n i)/a aims (f nl ) ~ 0.93 for the 
wavelet estimator. We have obtained a similar ratio for the 
bispectrum estimator in the same ideal experiment using the 
Fisher matrix and Gaussian simulations. 

Finally, we have studied the dependence of the wavelet 
estimator on the angular scale. In Fig. [7] we present the 
o'Fisher(fni) obtained with the wavelet estimator vs the 
Rmin and the Rmax- We have used the PCA method in 
order to minimize the influence of the errors in the in- 



© 0000 RAS, MNRAS 000, 000-000 



On the optimality of the spherical Mexican hat wavelet estimator for the primordial non-Gaussianity 7 



jj 10 




500 1000 1500 2000 

Number of components 



20 
15 
10 

5 




$ 


1 


1 ' 


;| 


o-(f nl ) with the 


<r(f nl ) with wcvelets O 
bispectrum (l mn , = 1535) ■ . 


" O 




- 




\ 



500 1000 1500 2000 
number of components 



Figure 6. Left: list of sorted eigenvalues of the covariance matrix C among the list of 1771 third order statistics. Right: cr(fni) obtained 
with the wavelet estimator using different subsets of eigenvalues of the C matrix following the PCA method. The vertical line shows the 
limit where the ratio of the maximum and minimum considered eigenvalues is 10 12 . 
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Figure 7. Left: <r(/ n ;) vs the maximum angular scale R max considered for the analysis. Right: cr(f n i) vs the maximum angular scale 
Rmin considered for the analysis. 



verse of the C matrix by using only the eigenvalues with 



Di/Dr, 



< 10" 



We can see that (TFisher(fni) is highly de- 



pendent on Rmin ■ This dependence is explained because the 
small scales map better the higher multipoles (even for scales 
smaller than the pixel size), while the large scales are all cen- 
tered around low multipoles. We have also found that the 
most significant contribution to the local non-Gaussianity 
is given by combinations of small and large scales, which 
is in agreement with the squeezed configurations of the 
shape form of this kind of non-Gaussianity (see for e xample 
iFergusson fc Shellardll2009l ; Ijeong fc KomatsulbOOfjl ) . 



3.2 Estimated error bars for / n j with a 
WMAP-like experiment 

We have estimated the expected dispersion of the local f n i 
parameter for the WMAP 5-year data and WMAP 7-year 
data. In this analysis we have considered the instrumental 
properties of the experiment, the noise level and the recom- 
mended sky cuts of each data release. We have compared our 
results with the values obtained wi th the optimal estimator 
for WMAP 5-yr (|Smith et alj|2009l '). where a(fni) = 21, and 
with the optimal estimator for WMAP 7-yr (jKomatsu et al.l 
l2010l ) where a(f n i) = 21. 

For WMAP 5-yr analysis we have used 10 4 Gaussian 
simulations of the combined V+W WMAP 5-yr data to es- 



timate the covariance matrix C and 300 non-Gaussian sim- 
ulations to estimate the expected values of the third order 
moments. We have normalised these non-Gaussian simula- 
tions to the ACDM model that best fits the WMAP 5-year 
data. As we are using a fraction of the full sky, we need to ex- 
tend the masks in order to avoid the propagatio n of the zeros 
that a re masking the Galaxy to other regions (|Curto et al.l 
l2009aT ). specially for large angular scales. This imposes a 
limit on the largest scale available for the analysis. In par- 
ticular, from the set of scales given in Fig. [51 we use all 
the scales with at least 10% of the sky. Following the c rite- 
ria to extend the masks defined in lCurto et al.l ( 2009aT l . we 
have 15 available scales. Their corresponding masked area 
is given in Table [l] In Fig. [8] we plot the eigenvalues of the 
C matrix for this case and the cr(fni) for different subsets 
of eigenvalues of C using the PCA method. Using 10 12 as a 
safe limit for the ratio of the maximum and minimum eigen- 
values, we obtain a(f n i) = 22 for this data. Compared with 
the optimal estimator based on the bispectrum (|Smith et al.l 
120091 ). we obtain ver y similar results. C ompa red with the 
results presented by ICurto et alj (|2009al ) and ICurto et al.l 
l|2009q ). we obtain significant improvements. The better re- 
sults obtained here are explained because of the wider ratio 
of large-to-small angular scales considered now. This allows 
to have more third order moments with squeezed triangles 
for their corresponding three angular scales. These are the 
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Table 1. Considered angular scales and the fraction of the sky that is masked out by each extended mask for the WMAP 5-yr initial 
mask. Similar values are obtained for the WMAP 7-yr masks. 
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Figure 8. Left: list of sorted eigenvalues of the covariance matrix C among the list of 680 third order statistics (15 scales) for WMAP 
5-year data. Right: cr(/ n ; ) obtained with the wavelet estimator using different subsets of eigenvalues of the C matrix following the PCA 
method. The vertical line shows the limit where the ratio of the maximum and minimum considered eigenvalues is 10 12 . 



combinations where the most sig nificant part of the loca l 
non-Gaussianity signal is located (|Yadav fc WandeltJuOlfJ) . 
For the case of WMAP 7-year data, we have used 10 4 
Gaussian simulations of the combined V+W WMAP 7-yr 
data to estimate the covariance matrix C and the same 300 
non-Gaussian simulations normalised to the ACDM model 
that best fits WMAP 7-year data to estimate the expected 
values of the third order moments. We consider the same 15 
angular scales as for the WMAP 5-year case masking a simi- 
lar area. In Fig. [5] we plot the eigenvalues of the C matrix for 
this case and the cr(fnl) for different subsets of eigenvalues of 
C using the PCA method. Using 10 12 as a safe limit for the 
ratio of the maximum and minimum eigenvalues, we obtain 
<r(f n i) = 21 for this data map. Com pared with the optima l 
estimator based on the bispectrum (jKomatsu et alj|2010l ). 
we obtain equally optimal error bars with the SMHW. 



4 CONCLUSIONS 

We have developed an efficient method to constrain the lo- 
cal f n i with the CMB anisotropies based on wavelets. We 
have found the dependence of the third order moments de- 
fined in Eq.Uon /„; and the cosmological model through the 
primordial bispectrum (see Eg. I12[) . On the other hand we 
have found an analytical expression for the covariance ma- 
trix for all the third order statistics (see Eg.ll6[). Assuming 
a Gaussian-like distribution for the third order moments, 
we have estimated the variance of f n i through the Fisher 
matrix (see Subsect. 12. 4p . This variance is compared with 
the variance obtained with the same method for the bispec- 
trum in Sect. [3] Both cases have been applied to an ideal 



experiment with an angular resolution of 6.9 arcmin and 
without instrumental noise. After applying Principal Com- 
ponent Analysis in order to minimize the influence of the 
errors of the inversion of the covariance matrix, we have 
found that cr(fnl) = 5.4 when we use the wavelets while 
&(fni) = 5.4 for the bispectrum up to £ ma x = 1535. In ad- 
dition, considering the case of a more realistic experiment 
with anisotropic and incomplete sky, such as the WMAP 
data, we have obtained o{f n i) - 22 for V+W WMAP 5-year 
data and a(f„i) = 21 for V+W WMAP 7-year data. These 
results are almost equal to the values obtained with the op- 
tima l bispectrum estimator where a(f n i) = 21 (jSmith et al.l 
2009; Komats u et al.ll2010l ). All these results indicate that 
wavelets can be as efficient as the bispectrum to detect the 
non-Gaussianity of the local type. Apart from the efficiency 
of the tool (about 7 times faster than the bispectrum estima- 
tor for a WMAP-like experiment) it is remarkable that, as 
they are different statistical estimators, wavelets may be sen- 
sitive to different systematics in real data. Moreover wavelets 
allow us to test the i sotropic character of the /„; parameter 
l|Curto et all l2009bl ; iRudiord et all |2010| ) by studying dif- 
ferent regions of the sky. We stress the importance of this 
statistical tool as an efficient alternative to measure local 
f n i in experiments such as Planck. In forthcoming works we 
will apply this tool on WMAP real data to constrain the 
local and other configurations of the non-linear coupling pa- 
rameter /„;. 
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Figure 9. Left: list of sorted eigenvalues of the covariance matrix C among the list of 680 third order statistics (15 scales) for WMAP 
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